DAS utilizing GPU

Created by Paul Klippel in May 2022

---

This notebook will run and then plot using Delay and Sum photoacoustic image reconstruction method.

The Delay and Sum algorithm works by synchronizing the arrival times for an array of transducer elements, summing their values, and averaging. The result represents the initial pressure value for a discrete point in space. Repeated over a region in space and one can create an initial pressure image from the photoacoustic energy transfer function.

---

Import Packages

from scipy.io import loadmat
import numpy as np
import plotly.express as plt
import time as time
import math as mt
from numba import cuda

Kernel Function: InitIndexMat
Inputs:

xImageLocs: all x-values for each reconstructed pixel
yImageLocs: all y-values for each reconstructed pixel
sensorLocs: [x, y] position for each element
indexMat: empty matrix for each the corresponding index of each element's pixel
signalLength: length of rfdata's first dimension
sos: medium speed of sound
dt: inverse of the sampling frequency for rfdata

Output:

N/A

Purpose: To calculate the index matrix for each pixel and element. It assumes that indexMat will be used to index a flattened signal matrix

@cuda.jit
def InitIndexMat(xImageLocs, yImageLocs, sensXloc, sensYloc, indexMat, signalLength, conversion):
    layer, locx, locy = cuda.grid(3)
    if layer < indexMat.shape[0] and locx < indexMat.shape[1] and locy < indexMat.shape[2]:
        indexMat[layer, locx, locy] = (mt.sqrt((xImageLocs[locx, locy]-sensXloc[0,layer])**2+
                                              (yImageLocs[locx, locy]-sensYloc[0,layer])**2)/conversion
                                          + (signalLength*layer))
    cuda.syncthreads()

Class Name: VariableSetup
Inputs:

filepath: exact file location for the .mat datafile
xlims: min and max values for the x-axis in the reconstruction region  
ylims: min and max values for the y-axis in the reconstruction region  
res: pixel resolution in the reconstruction region  

Attributes(memory host):

.sens_count(CPU): number of transducer elements  
.rfdataALL(CPU): If more than one frame is present this holds all frames
.rfdata(GPU): time vector with pressure values for the first frame recorded by each element  
.sensLocs(CPU): [x, y] position for each element   
.c0(CPU): medium speed of sound  
.dt(CPU): inverse of the sampling frequency  
.xMat(CPU): all x-values for each reconstructed pixel  
.yMat(CPU): all y-values for each reconstructed pixel  
.allIndex(GPU): indeces corresponding to each pixel in the layered delay matrix for all elements  
    l x m x n, l - corresponding to each element layer, m x n - corresponding to the reconstruction region

Methods:

__init__(self, filepath, xlims, ylims, res):
    initialize a VariableSetup object with attributes above and one conditional;
    For multiple frame rfdata, only the first frame is sent to the GPU.

SetupGPU(self, tb_shape):
    create the index matrix for all image reconstruction frames, calls above "InitIndexMat" kernel function on cuda.
    Inputs:
        self - VariableSetup object
        tb_shape - threads per block shape for cuda kernels
    Return:
        blocks_grid - shape of the blocks per grid for future image reconstruction

Purpose: To hold the necessary information for calculating DAS images given a constant system setup and speed of sound.

class VariableSetup:
    def __init__(self, filepath, xlims, ylims, pixLength):
        dat = loadmat(filepath)
        self.sens_count = np.squeeze(dat['ele'].astype(int))
        
        if np.array(dat['rfdata']).ndim > 2:
            self.rfdataALL = np.array(dat['rfdata'])
            self.rfdata = cuda.to_device(self.rfdataALL[:,:,0])
        else:
            self.rfdata = cuda.to_device(np.array(dat['rfdata']))
            
        self.sensXlocs = np.array(dat['x0']) * 1e-3
        self.sensYlocs = np.array(dat['z0']) * 1e-3
        self.c0 = dat['sos'] * 1e3
        self.dt = 1 / (dat['fs'] * 1e6)
        self.xMat = np.transpose(np.tile(np.linspace(xlims[0], xlims[1], pixLength), (pixLength, 1)))
        self.yMat = np.tile(np.linspace(ylims[0], ylims[1], pixLength), (pixLength, 1))
        self.allIndex = cuda.device_array((int(self.sens_count),pixLength,pixLength),dtype = 'uint64')
        
    def SetupGPU(self, tb_shape):
        threads_block = tb_shape
        blocks_gridx = mt.ceil(self.allIndex.shape[0]/threads_block[0])
        blocks_gridy = mt.ceil(self.allIndex.shape[1]/threads_block[1])
        blocks_gridz = mt.ceil(self.allIndex.shape[2]/threads_block[2])
        blocks_grid = ( blocks_gridx, blocks_gridy,blocks_gridz)
        InitIndexMat[blocks_grid, threads_block](self.xMat,self.yMat,self.sensXlocs,self.sensYlocs,
                                                 self.allIndex,int(self.rfdata.shape[0]),float(self.c0*self.dt))
        return blocks_grid

Kernel Function: FrameImaging
Inputs:

indexArr: Array from VariableSetup class .allIndex
valueArr: Flattened array from VariableSetup class .rfdata, only one "frame" at a time.
outputArr: Array with the same shape and size as the indexArr, used to store initial pressure values for each pixel

Output:

N/A

Purpose: To calculate the initial pressure values for each pixel averaged by each transducer element.

@cuda.jit
def FrameImaging(indexArr, valueArr, ImgArr):
    layer, locx, locy = cuda.grid(3)
    if layer < indexArr.shape[0] and locx < indexArr.shape[1] and locy < indexArr.shape[2]:
        ImgArr[layer, locx, locy] = valueArr[indexArr[layer, locx, locy]]/indexArr.shape[0]
    cuda.syncthreads()

Main Function

---

Setup Variables- This is only done for the first frame, setting up the index matrix

xlimits = [-0.015, 0.015]
ylimits = [-0.015, 0.015]
pixNum = 960
file = "D:\Documents\PSU\Research\DelayAndSum\Circ_Array\Leaf_exp\leaf.mat"

DASvars = VariableSetup(file,xlimits,ylimits,pixNum)
tpb = (4,16,16) #Unique to GPU, mine fits this shape well
bpg = DASvars.SetupGPU(tpb)

imageGPU = cuda.device_array((int(DASvars.sens_count),pixNum,pixNum), dtype = float)

C:\Users\paul\anaconda3\lib\site-packages\numba\cuda\cudadrv\devicearray.py:790: NumbaPerformanceWarning: Host array used in CUDA kernel will incur copy overhead to/from device.
  warn(NumbaPerformanceWarning(msg))

Create the image- This can be run in a loop with multiple frame rfdata, re-assigning rfdata for each frame

GPUstart = time.time()
FrameImaging[bpg, tpb](DASvars.allIndex, DASvars.rfdata.ravel(order='F'), imageGPU)
onGPUend = time.time()
image = imageGPU.copy_to_host()
offGPUend = time.time()
print('Time taken on GPU: ', onGPUend-GPUstart, ' [s]')
print('Time taken total:  ', offGPUend-GPUstart, ' [s]')

Time taken on GPU:  0.20968008041381836  [s]
Time taken total:   1.0188746452331543  [s]

Show the Reconstruction- optional to view the reconstructed image

fig = plt.imshow(np.mean(image,axis=0), 
                 labels=dict(x='mm', y='mm'), x=(DASvars.xMat[:, 0]).astype(str), y=(DASvars.yMat[0, :]).astype(str))
fig.show()

---

END OF MAIN PROGRAM

---

Below are code blocks for time comparison between CPU and GPU implementations

GPU based function

This section records the time required for 10 iterations of the gpu accelerated reconstruction.

timesNoCopy = np.zeros(10)
timesCopy = np.zeros(10)
for i in range(10):
    start = time.time()
    FrameImaging[bpg, tpb](DASvars.allIndex, DASvars.rfdata.ravel(order='F'), imageGPU)
    timesNoCopy[i] = time.time()-start
    image = imageGPU.copy_to_host()
    timesCopy[i] = time.time()-start

CPU based function

This section records the time required for 10 iterations of the cpu based reconstruction.

Reconstruction agreement has already been verified elsewhere.

class VariableSetup:
    def __init__(self, filepath, xlims, ylims, res):
        dat = loadmat(filepath)
        self.sens_count = np.squeeze(dat['ele'].astype(int))
        self.rfdata = np.array(dat['rfdata'])
        self.sens_xlocs = np.array(dat['x0']) * 1e-3
        self.sens_ylocs = np.array(dat['z0']) * 1e-3
        self.c0 = dat['sos'] * 1e3
        self.dt = 1 / (dat['fs'] * 1e6)
        xlength = int((xlims[1]-xlims[0])/res+0.5)
        ylength = int((ylims[1]-ylims[0])/res+0.5)
        self.xMat = np.transpose(np.tile(np.arange(xlims[0], xlims[1], res), (ylength, 1)))
        self.yMat = np.tile(np.arange(ylims[0], ylims[1], res), (xlength, 1))
        self.allIndex = np.zeros((self.sens_count, xlength, ylength), dtype='uint32')

    def initialize_index(self):
        shift = 0
        for i in np.arange(0, self.sens_count):
            self.allIndex[i, :, :] = np.round(np.sqrt(np.square(self.xMat - self.sens_xlocs[0, i]) +
                                         np.square(self.yMat - self.sens_ylocs[0, i])) / self.c0 /
                                              self.dt)+(i*len(self.rfdata))
            shift += len(self.rfdata)


'''
Main
'''
xlimits = [-0.015, 0.015]
ylimits = [-0.015, 0.015]
dx = 31e-6

setup_vars = VariableSetup("D:\Documents\PSU\Research\DelayAndSum\Circ_Array\Leaf_exp\leaf.mat", xlimits, ylimits, dx)
setup_vars.initialize_index()

timesCPU = np.zeros(10)
for i in range(10):
    start = time.time()
    flat_rfdata = setup_vars.rfdata.flatten('F')  # flattened rfdata for linear indexing
    img = np.mean(flat_rfdata[setup_vars.allIndex.astype(int)], 0)
    timesCPU[i] = time.time()-start

print("Average on GPU w/o copy: %g [s]" % np.mean(timesNoCopy))
print("Average on GPU w/  copy: %g [s]" % np.mean(timesCopy))
print("Average time on CPU:     %g [s] \n" % np.mean(timesCPU))
print("Fastest on GPU w/o copy: %g [s]" % np.min(timesNoCopy))
print("Fastest on GPU w/  copy: %g [s]" % np.min(timesCopy))
print("Fastest time on CPU:     %g [s]" % np.min(timesCPU))
print("-------------------------------------------------------")
factorWO = np.mean(timesCPU)/np.mean(timesNoCopy)
print("Time Factor w/o copying: %g " % factorWO)
factorW = np.mean(timesCPU)/np.mean(timesCopy)
print("Time Factor w/  copying: %g \n" % factorW)
diffWO = np.min(timesCPU)-np.min(timesNoCopy)
print("Fastest time diff w/o:   %g [s]" % diffWO)
diffW = np.min(timesCPU)-np.min(timesCopy)
print("Fastest time diff w/ :   %g [s]" % diffW)

Average on GPU w/o copy: 0.00039947 [s]
Average on GPU w/  copy: 0.966129 [s]
Average time on CPU:     2.84349 [s] 

Fastest on GPU w/o copy: 0 [s]
Fastest on GPU w/  copy: 0.886762 [s]
Fastest time on CPU:     2.77889 [s]
-------------------------------------------------------
Time Factor w/o copying: 7118.16 
Time Factor w/  copying: 2.94318 

Fastest time diff w/o:   2.77889 [s]
Fastest time diff w/ :   1.89213 [s]

Future improvements:

• Setup shared memory for averaging the output image on the GPU
• Setup plotting directly from GPU for stream optimization
• Create frame buffer system for realtime imaging